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Abstract 

The quantum-classical Liouville equation describes the dynamics of a quantum subsystem coupled to 
a classical environment. It has been simulated using various methods, notably, surface-hopping schemes. 
A representation of this equation in the mapping Hamiltonian basis for the quantum subsystem is derived. 
The resulting equation of motion, in conjunction with expressions for quantum expectation values in the 
mapping basis, provide another route to the computation of the nonadiabatic dynamics of observables that 
does not involve surface-hopping dynamics. The quantum-classical Liouville equation is exact for the spin- 
boson system. This well-known model is simulated using an approximation to the evolution equation in the 
mapping basis and close agreement with exact quantum results is found. 
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I. INTRODUCTION 



Nonadiabatic quantum mechanical effects are known to be important for the description of the 
dynamics of many chemical and biological processes. Photochemical dynamics, proton and elec- 
tron transfer reactions and vibrational relaxation processes are just a few examples where quantum 
effects play significant roles. Due to the difficulty of simulating the full quantum dynamics of 
large, complex, many-body systems, various mixed quantum-classical and semiclassical schemes 
have been developed. Here we consider quantum dynamics based on the quantum-classical Liou- 
ville equation (see Ref. yj] and references therein). 



^Pw{R,P,t) = -j[Hw,Pw{t)] (1) 
at n 

+^i{Hw.Pwit)}-{pwit).Hw}), 

where [A,B\ is the commutator and is the Poisson bracket for any operators A and B. The 

density matrix pwiR^Pj) is a function of the environmental phase space variables {R,P) and is 
an operator in the degrees of freedom of the quantum subsystem. The Hamiltonian Hw includes 
terms describing the quantum subsystem, its environment and the coupling between these parts of 
the system. This equation has been used to the describe nonadiabatic dynamics on coupled elec- 
tronic states^, vibrational dephasing'^, proton transfer reactions^'^'^i'^'^ and population relaxation 
in the spin-boson model^^'^^ to name a few examples. The simulation of the dynamics using this 
equation presents challenges and a number of different schemes have been devised for this pur- 
pose. Often the simulation methods are based on specific representations of the quantum degrees 
of freedom. For example, surface-hopping dynamics that make use of the adiabatic basis have 
been constructed^^*^*^, evolution of the density matrix in the diabatic basis has been carried out 
using a trajectory-based algorithm^i and a representation of the dynamics in the force basis has 
been simulated using the multithreads algorithmic*^. 

The discrete quantum degrees of freedom of the system can be described by the "classical elec- 
tron analog" model^ or the mapping formalis m^^'^^ i^t^. Extending Schwinger's angular momen- 
tum formalism^*^ to the A^-level case, the mapping formulation employs a quantum-mechanically 
exact mapping of discrete electronic states onto continuous variables; thus, the dynamics of both 
electronic and nuclear degrees of freedom are described by continuous variables^. The mapping 
basis has been used to compute quantum dynamics in the context of semiclassical path integral 
formulations of the theor5*i^ii^ii^i^ and in linearized path integral methods22j21)2^. in this paper 



we show how the quantum-classical Liouville equation can be written in this mapping representa- 
tion. The resulting evolution equation, like the basis-free quantum-classical Liouville equation ([1]) 
from which it was derived, provides a useful description of the dynamics of a quantum subsystem 
coupled to its environment. Since the quantum-classical Liouville equation is exact for any quan- 
tum system bilinearly coupled to a harmonic bath, so is its representation in the mapping basis 
presented here. The spin-boson model is of this type and this standard test model, for which exact 
quantum results are available, is employed to illustrate features of the simulation of the mapping 
form of the quantum-classical Liouville equation. In particular, we show that an approximation to 
the evolution operator allows one to accurately simulate the evolution using few trajectories with 
an algorithm that does not involve surface hopping dynamics. Comparisons with the results of 
other simulation algorithms are made. A discussion of the applicability of this representation of 
the theory to general many-body quantum systems is given in the last section of the paper. 



n. QUANTUM-CLASSICAL DYNAMICS IN THE MAPPING BASIS 

We consider a quantum mechanical system that is partitioned into a subsystem and bath. The 
expectation value of an operator B{t) can be written generally as 

B{t) =Tr B{t)p =Tr' J dX BwiXj)pwiX), (2) 
where a partial Wigner transform over bath degrees of freedom, 

Bw{X) = I dZe''-^/\R-^\B\R + ^), (3) 

has been taken. Here, X = {R,P) denotes phase space variables of the bath. The initial density 
matrix is pw{X). The partially Wigner transformed Hamiltonian of the system can be written as 

^w^^ + ^+v,{4)+VB{R)+Vc-{q,R). (4) 

where the subscripts s, B and c denote the subsystem, bath and coupling, respectively. Letting 
hs = p^/2m + Vs{4) be the subsystem Hamiltonian, whose eigenvalue problem is = £^1-^), 
we can write the expectation value of B{t) in the form 

5(0 = £ / dXB^^'{Xj)pi;^{X). (5) 

in the subsystem basis. 



A. Average value in mapping basis 



Next, we write this expectation value in the mapping basis by noting that any opera- 
tor Bw{X) can be decomposed as B\y{X) = B^^' {X)\X.) The evolution of the A'^- 
state subsystem can be conveniently replaced, using mapping relations, with that of N ficti- 
tious harmonic oscillators with occupation numbers limited to or 1, namely, |A) \mx) = 
|Oi,---,l;i,,- ■ ■ 0„)^^idMldM2iiii^iL2i^i2dail, The matrix element of an operator may then be written 
in the mapping form, B^^' (X) = {?i\Bw{X)\?i') = {m^\B,„{X)\mi,), where 

B^{x) = Y,B^^'{x)dldx,. (6) 

The mapping annihilation and creation operators are given by 




ax = \l i;r{qx+iPx): a\ = \l^_{qx-ipx), (7) 



and satisfy the commutation relation [dx , = 5xx' ■ Explicitly, we may write 

= i^mqx' +PXPX' -Apx^x' -qxPx')\- (8) 

One may easily verify that the matrix elements of B„-i{X) in the mapping basis are identical to 
those of Bw{X) in the subsystem basis. 

In the analysis that follows it is convenient to work in a Wigner representation of the mapping 
basis. To this end we introduce a coordinate representation of the mapping states and annihilation 
and creation operators, 

{q\mx) = (<?i,^2,---,^yv|0i,...,lA,---,0A^) 

= ^o{gi)-hiqx-i)M(ix)---^oM (9) 

and 

with an analogous expression for {q\d'^^\q'). Here 

^,{qx) = V2{7thY'^'qxe-^l/^\ (11) 
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Equation ^ may be written in the mapping basis using the coordinate representation to obtain 

5(7) =£ [dX{m^\B„{X,t)\m^,) 



X {tnx'\pmiX)\mx) 



= ^ dX dqdq'dq'dq" (mx\q){q\Bm{X,t)\q') 

{q'\m^,){m^,\q"){q"\p^{X)\q"'){q"'\m^). (12) 

Note that the coordinate space dimension of the mapping variables is N. We may now introducing 
the Wigner transforms of the coordinate space matrix elements of the mapping variables, 



(r-||5„(X,0k + |) = / dpe-'P^^/'^B,^{x,X,t) 

(/-||p„(X)|/+|) = j dp'e-^P'^'/'^pUx^X), (13) 



where x = {r,p) are the phase space coordinates of the mapping variables. Using these definitions 
Eq. (fT2)) can be written as 

B{t) = jdXdxdxB,n{x,X,t)f{x,x)pm{x,X) 

= j dXdxB,n{x,X,t)p„^{x,X), (14) 

where pm{x,X) = J dx' f{x,x')pm{x' ,X) and 



X - |)(/ + ^|m;,).-'(^-+^'-')/l (15) 



The function f{x,x') can be computed explicitly using Eqs. dH) and (fTTjl . 



B. Evolution equation in mapping basis 

In quantum-classical dynamics the time evolution of an operator may be described by the 
quantum-classical Liouville equation,— 

jBwit) = ^[HwJwit)] (16) 
-^{{Hw.Bwit)}- {Bw{t),Hw}). 
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In order to make use of Eq. (1141) we must cast this equation in the mapping basis. 

Using the results of the previous subsection, the matrix elements of any operator Cw can be 
written in the form, 

(A|Cm/|A') = j dqdq'{mx\q){q\C,n\q'){q'\mx'). (17) 
Furthermore, if Cw = AwByy is a composition of operators, we have 

(A|Civ|A') = {mx\Cm\mx,) = {mx\AinB,n\mxf). (18) 

Using Eqs. (fTTI) and (fTSi) . we can write the quantum-classical Liouville equation as 

^^{q\B^{t)\q') = '-{qWH^^Bmmq') (19) 

~{{q\{H^,B^{t)}\q') - {q\{B^{t),HM))- 

Taking the Wigner transform over the mapping coordinate space we obtain 

^B^{x,X,t) = ^//„sin(^)5„(0 (20) 
at n 2 

dH,n i^^m. dBmit) , P dB„^{t) 



dR ' 2 ' dP M dR ' 
where the negative of the Poisson bracket operator on the mapping phase space coordinates is 
defined as A,„ = Vp ■ — ■ V^. In writing this equation we have used the fact that the Wigner 
transform of a product of mapping operators is given by 

{A^{X)B,n{X))w ^A^{x,X)e^^"'/^'B,,{x,X). (21) 

The Wigner transform of a mapping variable is given by 

AUx.X) = ^£A^^'(7?)(r;,a'+PAPA'-^5;,A'), (22) 
In particular the mapping Hamiltonian takes the form 

H^ix.X) = :^ + Vb{R) (23) 

XX' 

where hxx'{R) = {^\p^ /'2m + Vs{q) + Vc{q,R)\X') and we have used the fact that = h^'x- 
Given this form of the Hamiltonian one may show that 

H^KrnBm = T F h^x' [px 3^ " rx^—)B,n, (24) 



fi dri, dri dpi, dpx 

and 

HmAlBm = 0. (whenn > 3) (26) 



Then, using these relations, we can simplify Eq. (1201) to derive the quantum-classical Liouville 
equation in the mapping basis: 

/ P d dH,„ d \ 

h ^ dh^it d d d d d 

Since the quantum-classical Liouville equation is exact for an arbitrary quantum subsystem bilin- 
early coupled to a harmonic bath, the mapping version of this equation, Eq. (|27l) . is also exact for 
such systems. 

The first term in Eq. (|271) is the quantum evolution of the subsystem in the mapping phase 
space, while the second term describes the evolution of the bath where the forces involve the 
mapping coordinates. The complicated third term represents the higher-order correlations between 
the subsystem and the bath. The evolution equation can be written more compactly as 

—Bm{x,X,t) = —{Hm-,Bm{t)}xX (28) 

^-^L^,^^^,d7;^j^,j^^^-dp''-^'^ 

where denotes a Poisson bracket in the full mapping-bath phase space of the sys- 

tem. The last line of this equation defines the quantum-classical Liouville operator in the mapping 
basis, 

/^„ = /^o^iX, (29) 

where 

/J5f^ = -{H,n,B,nUx. and (30) 

, ^ nydhxx^ _d d_ d d_ d_ 

'"^^ 8^^, dR "^drydri^ dpx,dpx^' dp- 
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The complex form of the makes simulation of the dynamics in the mapping basis difficult. If 
the evolution equation is approximated by z^^, simulation of the dynamics in terms of Newtonian 
trajectories is straightforward in view of Poisson bracket form of the resulting equation of motion. 
The validity of this approximation must be determined for specific applications. In the next section 
we apply this equation to the spin-boson model and show that accurate results can be obtained 
when the last term on the right side of this equation is neglected. 

Ill, SPIN-BOSON MODEL 

The spin-boson model is often used as a test case for quantum simulations of many-body sys- 
tems and we present the results of simulations of this model using the quantum-classical Liouville 
equation in the mapping basis. The spin-boson model describes a two-level system bilinearly 
coupled to a harmonic bath of Nb oscillators with masses Mj and frequencies (Oj. The system 
Hamiltonian is given by 



where Hq = Y.j (Pj /2Mj +Mj(ojE?- /2j and y{R) = -Y.cjRj. The energy gap of the isolated 
two-state system is 2hQ.. From Eqs. (|23l) and (|3TI) we can obtain 



Equation (|28l) is an exact evolution equation for the spin-boson model. Previous simulations of the 
quantum-classical Liouville equation in the adiabatic basis have been carried out using a Trotter- 
based scheme and were able to reproduce the exact results for a wide range of system parameters.— 
Consequently, the results presented here can be viewed as a test of the utility of the simulation 
schemes that use the mapping basis to represent quantum-classical Liouville dynamics. 

As in previous studies we assume that the initial density matrix is uncorrelated so that the 
subsystem is in the ground state and bath is in thermal equilibrium, namely. 




(31) 




(32) 




(33) 



where the Wigner distribution of the bath Pb(X) is given b5*=2a2a 




(34) 
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with u'j = UjCOthuj and uj = j5h(Oj/2. The subsystem initial density matrix in the Wigner- 
transformed mapping basis is 

p,„,W = {2nn)-\2n)-\rj+pl-n). (35) 

Using the results in Eq. (fT4l) . the time evolution of the population difference between the ground 
and excited states, given by the expectation value of Pauli matrix a^, can be written as 

= 1^1 dXa,'''it)p^'H0)PB{X) (36) 

XX' 



dXdx Ozm {x,X,t) psm {x) Pb{X), 



where psm{x) has the explicit form 



^■M = ^M + pi-^^)e^^'"'-''^^'- (37) 
The initial value of the Wigner-transformed mapping representation of is o^m{x) = (r^ + 



To compute <yv{t) we need to solve for o^mix^Xj) using Eq. (l28l) . This equation is difficult to 
solve because of the structure of the last term of the quantum-classical Liouville operator in the 
mapping basis. For the spin-boson model one may show by direct calculation for short times that 
the last term does not contribute until the fifth-order initial derivative of Ozm{x). This suggests 
that it may be possible to obtain a useful approximate solution by neglecting the last term in the 
evolution equation (|28l) so that 

j^O,m{x,Xj) ^ i^^,0„n{t)- (38) 

The dynamical variable 0:^,n{x,X, t) evolves by Newtonian equations of motion and admits a solu- 
tion in terms of characteristics. The corresponding set of ordinary differential equations is 



X 

dR{t) P{t) dP{t) dH, 



m 



dt M ' dt dR{t)' 

Using this result, we obtain the simple form for the expectation value. 



(39) 



= (^) / dxdX pB{X)e-^^'+p'y'^ (40) 
^iri+p]-l){n{tf + pdtf-r2{tf-p2{t)'). 



The linear coupling in the spin-boson model is characterized by an Ohmic spectral density, 
J{co) = n'£cj/{2Mj(Oj)5{co- cOj), where cj = {^hAcoMjY/'^cOj, coj = -(Ofln(l - jAco/cOc) and 
Aco = (Oc (^l -g-^'na^/^'^j /Nb with cOc the cut-off frequency and ^ the Kondo parameter.— We 
used Nb = 20 and CO,nax = 4ft)f. Dimensionless units with time scaled by COc are used in the calcu- 
lations below. The equations of motion were integrated using the velocity Verlet algorithm with 
time step At = 0.1. 

The expectation value o^{t) in Eq. (I40l) may be computed by sampling initial bath and mapping 
variables from Gaussian distributions, reweighting to account for the form of the initial density 
matrix, and computing Ozin{x{t)). We have also carried out the calculations using focused initial 
conditions^^iil where the state mapping variables are initially taken to be ri = 1, pi = 1, r2 = 
and P2 = when state 1 is initially occupied. Both sampling methods yield comparable results 
but focused initial conditions require about a factor of ten fewer trajectories to obtain converged 
results for this model. 

We tested our method for the parameters for which numerically exact results are available. 
Approximately 10^ trajectories were used to obtain the results in the figures. Comparable results 
can be obtained even with ten times fewer trajectories. In Fig. [T] the results are compared for 
weak system-bath coupling with ^ = 0.09. The adiabatic energy gap is chosen as Q. = 0.4. For 
high temperatures, the time-dependent population difference exhibits incoherent behavior as in 
Fig. [2 a). Our results, as well as those of other methods such as LAND-map and LSC-IVR, 
show excellent agreement with the numerically exact results^^ for the high-temperature, weak- 
coupling case. The reproduction of the coherent or oscillatory behavior at low temperatures shown 
in Fig.[TJb) is a more severe test, especially at long times. Our results predict the correct frequency 
of oscillations but the magnitude of the oscillations are somewhat smaller at long times. 

In Fig. [21 we plot (cr^(O) for ^ rather high friction constant, ^ = 2, at high temperature of 
/3 =0.25. One can see that the accuracy of our results does not change for strong system-bath 
coupling and is consistently better than other approaches. 

As a final test, in Fig. [3] we show {Oz{t)) for two friction constants, ^ = O.l and 0.5, for a 
relatively low temperature, /3 = 3. The LAND-map approach predicts the slow incoherent decay 
instead of oscillation around zero. This discrepancy was attributed to the linearization approxi- 
mation which underestimates the coherent dynamics. Our results again show reliable accuracy 
both for weak and strong coupling. Our results are compared with those using the semiclassical 
influence functional formalism with four time slices. Similar accuracy is obtained. 
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IV. CONCLUSION 



The representation of the quantum-classical Liouville equation in the mapping Hamiltonian 
basis provides another way to simulate nonadiabatic dynamics. The complicated form of i^ln in 
Eq. (|29l ) in this basis leads to difficulties in the construction of simulation algorithms. If this term 
is neglected and the evolution is approximated by iS^^ the dynamics may be computed easily using 
an ensemble of trajectories. This is an excellent approximation for the spin-boson model and leads 
to a simulation scheme for nonadiabatic dynamics that does not involve surface-hopping. The 
extent to which this approximation is applicable to more general systems remains to be determined. 
The work also suggests that it may be possible to construct simulation algorithms that use evolution 
under as a zeroth order scheme about which corrections can be computed. The utility of the 
mapping formulation of the quantum-classical Liouville equation for the computation of general 
correlation functions is also another topic that is worth pursuing. 
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1 Electronic population difference {0^(1)) as a function of t for two dimensionless 
parameter sets: Q. = 0.4, E, = 0.09, and /3 = 0.25 (a) or 12.5 (b). The solid points 
are exact results^, the dashed lines are the LAND-map results^ and the dotted 
lines are the LSC-IVR resuks^^ 

2 Electronic population difference (o^it)) as a function of t for two parameter sets: 

= 2, /3 = 0.25, and Q. = 0.8 (a) or 1.2 (b). The solid points are exact results^, 
the dashed lines are the LAND-map results^-, the dot-dashed lines are the TDSCF 
results^ and the dotted lines are the LSC-IVR results^ 

3 Electronic population difference (c7j;(?)) as a function of t for two parameter sets: 
Q. = 1/3, P = 3, and ^ = 0.1 (a) or 0.5 (b). The solid points are exact results^—, 
the dashed lines are the LAND-map results^, and the open squares are the results 
obtained using imaginary time path integral semiclassical influence functional for- 
malism with four time slices^ 
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Figure 1: Electronic population difference (cTj.(f)) as a function of t for two dimensionless parameter sets: 
Q. = 0.4, = 0.09, and j8 = 0.25 (a) or 12.5 (b). The solid points are exact results^'', the dashed lines are 
the LAND-map results^! and the dotted lines are the LSC-IVR results^-. 
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Figure 2: Electronic population difference (cJ^(f)) as a function of t for two parameter sets: = 2, j8 = 0.25, 
and n = 0. 8 (a) or 1 .2 (b). The solid points are exact results^, the dashed lines are the LAND-map results^l, 
the dot-dashed lines are the TDSCF results^ and the dotted lines are the LSC-IVR results^^. 
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Figure 3: Electronic population difference (cJj.(f)) as a function of t for two parameter sets: Q.= 1/3, 
j8 = 3, and (^=0.1 (a) or 0.5 (b). The solid points are exact results^, the dashed lines are the LAND-map 
results^, and the open squares are the results obtained using imaginary time path integral semiclassical 
influence functional formalism with four time sUces^. 
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